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Much research has been done on the numerical integration of the radial 
Schrodinger equation: 



where ^^^^ ' is the centrifugal potential, V{x) is the potential, E is the 

energy and W{x) = ^^^^P + V{x) is the effective potential. It is valid 
that lim Vix) = and therefore lim W{x) = 0. 

x— >oo a;— >oo 

Many problems in chemistry, physics, physical chemistry, chemical 
physics, electronics etc., are expressed by equation (1). 

In this paper we will study the case of E' > 0. We divide [0, cxd] into 
subintervals [oj, hi] so that W{x) is a constant with value Wi. After this 
the problem (1) can be expressed by the approximation 
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y'l = {y/ — E)yi, whose solution is 

yi(x) = Ai exp (^Vw - E + Bi exp (-\/w -Ex), (2) 
A, Bi G M. 

There has been an extended bibhography on the development and 
analysis of numerical methods for the efficient solution of the Schrodinger 
equation: see for example [1]-[129]. 



2. Basic theory 

2.1. Explicit Runge-Kutta methods 

An s-stage explicit Runge-Kutta method used for the computation of 
the approximation of yn+i{x), when yn{x) is known, can be expressed 
by the following relations: 



Vn+l = Vn + ^hh 
1=1 

I \ 

ki = hf [xn + Cih, yn + Qij kj \ , i = 1, . . . , s (3) 

where in this case / {x,y{x)) = {W{x) — E) y{x). 

Actually to solve the second order ODE (1) using first order numer- 
ical method (3), (1) becomes: 

z'{x) = {W{x) - E) y{x) 

y'{x) = z{x) ^ ' 

while wc use two sets of equations (3): one for yn+i and one for Zn+i- 
The method shown above can also be presented using the Butcher table 
below: 














0-21 






C3 


031 


032 




Cs 


asi 


as2 ■ 


■ (ls,s-l 




hi 


b2 . 


■ bs-1 bs 



(5) 



Coefficients C2, . . . , must satisfy the equations: 
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i-l 

Ci = ^aij, i = 2,...,s (6) 

DEFINITION 1. [3] A Runge-Kutta method has algebraic order p 
when the method's series expansion agrees with the Taylor series ex- 
pansion in the p first terms: y^'^\x) = y^p.{x), n = 1,2, ... ,p. 

A convenient way to obtain a certain algebraic order is to satisfy a 
number of equations derived from Tree Theory. These equations will 
be shown during the construction of the new methods. 

2.2. Phase-Lag Analysis of Runge-Kutta Methods 

The phase-lag analysis of Runge-Kutta methods is based on the test 
equation 

y' = luy, ujeR (7) 

Application of the Runge-Kutta method described in (3) to the scalar 
test equation (7) produces the numerical solution: 

yn+i = aJlyn, a* = As{v^) + ivBs{v^), (8) 

where v = uh and As, Bg are polynomials in f ^ completely defined by 
Runge-Kutta parameters ajj, hi and Cj, as shown in (5). 

DEFINITION 2. [1] In the explicit s-stage Runge-Kutta method, 
presented in (5), the quantities 

t{v) =v — arg[a*(v)]; a{v) = 1 — |a*(?;)| 

are respectively called the phase-lag or dispersion error and the dissi- 
pative error. If t{v) = 0{v'^~^^) and a{v) = 0{v^'^^) then the method is 
said to be of dispersive order q and dissipative order r. 



3. Construction of the new trigonometrically fitted 
Runge-Kutta methods 

We consider the explicit Runge-Kutta method with 3 stages and 3rd 
algebraic order given in table (9). 
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(9) 



We will construct two new optimized methods. 

3.1. First optimized method with zero phase-lag 

In order to develop the new optimized method, we set free 63, while all 

other coefficients arc borrowed from the classical method. We want the 
phase-lag of the method to be null, so we satisfy the equation PL = 0, 
while solving for 63, where 



PL = 1/6 (6 + (-2-6 63)w^)tan(w)+w'^63 + l/6 (-5-663)1; 

So 63 becomes 

^ —6 tan {v) + 2 tan {v) f ^ + 5 u 
Qv {v tan {v) — v"^ + 1) 

and its Taylor series expansion is 

1 1 4 4 . 17 o 206 in 7951 ,2 

63 = -\ H v^^ H v^^ - ... 

6 30 315 2835 31185 12162150 

where v = ujh, is a real number and indicates the dominant frequency 
of the problem and h is the step-length of integration. 

3.2. Second optimized method with zero phase-lag and 
derivative 

As for the development of the second optimized method, we set free 
62 and 63, while all other coefficients are borrowed from the classical 
method. We want the phase-lag and its first derivative of the method 
to be null, so we satisfy the equations {PL = 0, PL' = 0}, while solving 
for 62 and 63, where 

PL = 1/6 (6 + (-362 - 663) v^) tan (v) + t; (-1/6 - 62 - ^3 + M^) 
PL' = -V tan (v) 62 - 2 u tan (v) 63 + 5/6 + (tan {v) f - 1/2 ^^62 
-l/2v%2 (tan {v)f + 2 bsv'^ - h^v"^ (tan {v)f - h - h 

(10) 

Then we have 
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^ 1 12 i;+i''' + l.aii(r)r--12 t.aii(r) — i,'-''(la,ii(i;))" 

2 6 ^— 3 t)-(-tan(t))+D(tan(D))"^ — tan(i')i''^+t)'^+t)'* (tan(t)))'^) 

7 _ 1 5t;3(tan(D))^-|-7?;3-19 tan(t;)j;2+6t;(tan(?;))^-6t;+6 tan(ti) 

^ ~ 3 v'^[-3v+ta.n{v)+v{ta.n(v))'''-ta.n{v)v'^+v^+v^(ta.n{v)f) 

The Taylor series expansion of the coefHcients are given below: 



h — 2 _ ^ „,2 _ ^ ,,4 _ 3526 „6 _ 173788 _,8 _ 354768808 ,,10 



''2—3 15 315 14175 467775 " 638512875 

h — 1 _L A „2 I ^ „4 I 4201 6 I 207349 8 , 423287713 
■^3 6 15 '-^ 126 ^ 14175 467775 638512875 



(12) 

where v = uh, is a real number and indicates the dominant frequency 
of the problem and h is the step-length of integration. 



4. Algebraic order of the new methods 

The following 4 equations must be satisfied so that the new meth- 
ods maintain the third algebraic order of the corresponding classical 
method (9). The number of stages is symbolized by s, where s = 4. 
Then we are presenting the Taylor scries expansions of the remainders 
of these equations, that is the difference of the right part minus the left 
part. 

1st Alg. Order (1 equation) 3rd Alg. Order (4 equations) 

tbi = l t hi4 = \ 

i=l i=l 

s 

2nd Alg. Order (2 equations) X) haijCj = \ 

s ^'^^^ 
i=l 

(13) 

4.1. Equations remainders for the first method 

We are presenting Rem which is the remainder for all four equations 
for the first method: 

i?em = -^^4-3^.6 + _17_^8 + _206_„10 + ... (14) 

4.2. Equations remainders for the second method 

The four remainders of the equations for the second method are: 
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Rem, = J- + J_ I 113. ^8 , 7rn_ lo , 
nemi u -|- 21 ^ 1575 ^ 66825 ^ ^ ' • • 

7?pmo - J- „2 _L ^ „4 I 2438 6 , 24091 8 , 245903309 10 , 
/ie/(t2 — 15 ^ 630 ^ 14175 ^ ^ 93555 ^ 638512875 ^ ■ ■ ■ 

Rem-i - i- 7,2 + 11 ^,4 + 2213 6 , 54634 8 , 37177279 10 , 
ntirna — iq O ^ 70 ^ 9450 " ^ 155925 ^ 70945875 " ^ " " " 

Rpm . - A 7,2 I ^ „4 I 4201 6 , 207349 8 , 423287713 10 , 
/ie/At4 — 15 f 126 ^ 14175 " ^ 467775 " ^ 638512875 " ^ ■ ■ ■ 

(15) 

We see that the two optimized methods retain the third algebraic 
order, since the constant term of all the remainders is zero. 



5. Numerical results 
5.1. The inverse resonance problem 

The efficiency of the two new constructed methods will be measured 
through the integration of problem (1) with ^ = at the interval [0, 15] 
using the well known Woods-Saxon potential 



Uq 

uq = —50, a = 0.6, xq = 7 and ui = 

a 

and with boundary condition y(0) = 0. 

The potential V{x) decays more quickly than ^-^^-j^, so for large x 
(asymptotic region) the Schrodingcr equation (1) becomes 

y"ix) = (^-^ - y{x) (17) 

The last equation has two linearly independent solutions kxji{kx) 
and kxni{kx), where ji and ni are the spherical Bessel and Neumann 
functions. When x ^ 00 the solution takes the asymptotic form 

y{x) ^ Akxji{kx) - Bkxni{kx) , . 

D[sin{kx-Trl/2)+taji{Si) cos{kx-TTl/2)], ^ ^ 

where Si is called scattering phase shift and it is given by the following 
expression: 

, yjxi) Sjxj+i) - yjxj+i) Sjxj) 

tan(d;) = — ^ — — r -. r—. r-, (19) 

y{xi+i) C{xi) - y{xi) C{xi+i) 
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where S{x) = kxji{kx), C{x) = kxni{kx) and Xi < Xj+i and both 
belong to the asymptotic region. Given the energy we approximate the 
phase shift, the accurate value of which is 7r/2 for the above problem. 

We will use three different values for the energy: i) 989.701916, ii) 
341.495874 and iii) 163.215341. As for the frequency iv we will use the 
suggestion of Ixaru and Rizea [2] : 



5.2. Nonlinear Problem 

y" = -lOOy + sin(y), with y{0) = 0, y'(0) = 1, t G [0,20 7r], 
y(207r) = 3.92823991 • 10"^ and w = 10 as frequency of this problem. 

5.3. Comparison 

We present the accuracy of the tested methods expressed by the 
— logiQ (error at the end point) when comparing the phase shift to 
the actual value 7r/2 versus the log^g (total function evaluations). The 
function evaluations per step are equal to the number of stages of 
the method multiplied by two that is the dimension of the vector of 
the functions integrated for the Schrodinger {y{x) and z{x)). In Figure 
1 we use E = 989.701916, in Figure 2 E = 341.495874 and in Figure 3 



We compare the two optimized methods and the corresponding classical 
explicit Runge-Kutta method for the integration of the Schrodinger 
equation and the Nonlinear problem. We see that the second method 
with the phase-lag and its first derivative nullified is the most efficient 
in all cases, followed in terms of efficiency by the optimized method 
with zero phase-lag and then by the corresponding classical method. 
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6. Conclusions 
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Schrodinger Equation (E = 989.701916) 




Figure 1. Efficiency for the Schrodinger equation using E = 989.701916 
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